Dispersive spectrum and orbital order of spinless p-band fermions in an optical lattice 
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We study single-particle properties of a spinless p-band correlated fermionic gas in an optical 
lattice by means of a variational cluster approach (VGA). The single-particle spectral function is al- 
most flat at half-filling and develops a strongly dispersive behavior at lower fillings. The competition 
between different orbital orderings is studied as a function of filling. We observe that an "antiferro- 
magnetic" orbital order develops at half-filling and is destroyed by doping the system evolving into 
a disordered orbital state. At low filling limit, we discuss the possibility of "ferromagnetic" orbital 
order by complementing the VGA result with observations based on a trial wave function. We also 
study the behavior of the momentum distribution for different values of the on-site interaction. 
Finally, we introduce an integration contour in the complex plane which allows to efficiently carry 
out Matsubara-frequency sums. 

PACS numbers: 03.75.Ss,71.10.Fd,79.60.-i,03.75.Lm 



I. INTRODUCTION 

Ultracold atomic gases in optical lattices constitute a 
promising system to simulate and investigate strongly 
correlated quantum phases as a function of their model 
parameters, which can be controlled experimentally in 
a large range This field of research has greatly ex- 
panded after the pioneering realization of the Superfluid 
to Mott-insulator transition by loading bosonic atoms to 
the lowest band of optical lattices 0, Q . Recent progress 
in this field has been achieved by loading the atoms in 
the first excited band, which makes the study of orbital 
physics possible in these systems [1, HI, Q • Orbital degrees 
of freedom play an important role in many solid-state 
materials: Many interesting phenomena such as metal- 
insulator transitions Q, superconductivity H | , colossal 
magnetoresistance @, half metallicity [13 EL ^^c, are 
rooted in the coupling of orbital with the other degrees 
of freedom (spin, charge, and phonon). The study of 
orbital physics in optical lattices, in a pure and tunable 
environment, is believed to be of great help to understand 
the complicated orbital issue of solid-state systems. 

The basic physics of cold atoms in the first ex- 
cited band can be captured by a p-band Hubbard like 
Hamiltonian 0]. Many novel phenomena and quan- 
tum phases have been predicted for the p-band bosons 
d, i, [13, [H, [ii, [3, e^, quantum stripe order [l3|, 
Wigner crystallization and bond algebraic liquid 

phase [l^. For the spin-1/2 p-band fermions, an anti- 
ferromagnetic order was found at half filling in both the 
strong and weak interaction regimes [l6j . and a robust 
ferromagnetic order was shown to exist for a large range 
of interaction and at band filling lower than half-filling 
p7| . Itinerant fcrromagnetism was also proposed in the 
honeycomb lattices in Ref. [ll] . Experimentally, the poj 
ulation of higher band was studied by Browaeys et al. 
and Kohl et al. for bosons and fermions respectively. 
Recent experiments performed by Miiller et al. were able 
to realize long lifetime p-band orbital bosonic systems [1] . 

In particular, the orbital exchange physics in the Mott 



state of an orbital-only model, which is realized by load- 
ing the single-component (spinless) fermions into p-band 
optical lattices (see the Hamiltonian in Eq. ([1])), has been 
studied by Zhao et al. [2l| and Wu 0| for various ge- 
ometry lattices. In these works, a new orbital exchange 
mechanism was found, and long-range orbital order was 
predicted. At the same time, a similar orbital-only model 
was proposed to describe the ferromagnetic plane in tran- 
sition metal oxides with to„ orbital degeneracy, such as 
Sr2VOi and K2CUF4 [23|. The spectral properties of 
this model in the half-filled case have been studied, and 
it was shown that a hole in a background of antiferromag- 
netic orbital order does not localize but moves coherently 
due to an effective three-site hopping term. 

Motivated by these previous works, we study numeri- 
cally this spinless p-band model on a square lattice with 
an emphasis on the excitation spectrum and orbital order 
away from half-filling. The paper is organized as follows: 
In Sec. [U we present the Hamiltonian of the model, and 
wc briefly summarize the method used to approximately 
solve it, namely, the variational cluster approach. As a 
byproduct, in this work, we introduce and adopt a more 
efficient method to carry out sums over Matsubara fre- 
quencies, which could also be applied to other problems. 
Details are given in Appendix \X[ In Sec. IIIIl we present 
the calculated results including the single-particle spec- 
trum, orbital orders, and momentum distribution. Fi- 
nally, wc draw our conclusions in Sec. IIVI 



II. HAMILTONIAN AND METHOD 

Wc consider an anisotropic 3D optical lattice with op- 
tical trapping frequency ujz ^ oJx — '-^y, so that the dy- 
namics in the z direction is essentially suppressed. Sup- 
posing that the lowest s orbital of the optical lattice is 
fully occupied by fermions, the other particles can only 
fill the degenerate Px and Py orbitals |21|, i23l . A fermionic 
gas, which is polarized into a single hyperfine spin state 
by magnetic field and loaded in such optical lattice, can 
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potential between the cluster and the lattice. The "op- 
timal" value for these variational parameters is deter- 
mined, in the framework of Self-energy Functional Ap- 
proach (SFA) [Ml, by requiring that the SPA grand- 
canonical potential 



FIG. 1: (Color online) Reference cluster for the VGA calcu- 
lation (a) consisting of L = 10 sites (gray). Schematic repre- 
sentation of "antiferromagnetic" (b), and of "ferromagnetic" 
orbital orders (c). 



be described by the following 2D spinless p-band Hub- 
bard Hamiltonian 

H = tap (4,aCr+e^,a + H .C^ + XI "r,xrir,2,(l) 

r a.J3—x.y r 

Here, ^ creates a fermionic atom in the pa orbital 
at position r, is the unit vector of the (3 direc- 
tion (the lattice spacing is set equal to unity), tap ~ 
t\\Sa[3 + t±{l — Sap) is the hopping amplitude of orbital 
Pa in direction (3, and U is the on-site repulsive inter- 
action between atoms in different orbitals. The longitu- 
dinal hopping i|| is positive but the transverse hopping 
t± is negative because of the odd parity of p orbitals 
dll, 113. In general, one has \t\\ \ 3> \t±\ for the strongly 
anisotropic shape of p orbital, therefore, we choose a typ- 
ical small value t± ~ — 0.05t|| in our calculation [3, 
and set i|| as the energy scale. Obviously, in the limit 
— > the hopping will be restricted to one dimension 
and the number of particles in a given orbital (pa) will 
be conserved for each chain oriented along a. There is 
no s-wave scattering for atoms in a single hyperfine spin 
state because of the Pauli exclusion principle. Therefore, 
the interaction U between atoms mainly comes from p- 
wave scattering, whose strength can be tuned using a 
p-wave Peshbach resonance [21| . It is argued that U can 
be increased to the order of the recoil energy En in the 
present experiment [23 |. 

The Variational Cluster Approach (VGA) [H, ^ 
is an extension of Cluster Perturbation Theory (CPT) 
[13, m, [H]. Within CPT, the original lattice is di- 
vided into a set of disconnected clusters, and the inter- 
cluster hopping parameters are treated perturbatively. 
Within VCA. additional ("virtual") single-particle terms 
are added to the cluster Hamiltonian, to obtain a so- 
called reference system, and then subtracted perturba- 
tively. (So that if the perturbative treatment was exact, 
results would not depend on these terms). These single- 
particle terms can contain "Weiss" fields to describe a 
particular ordered state, but also other Hamiltonian pa- 
rameters, such as, for example, an offset in the chemical 
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is stationary within this set of variational parameters. 
Here, Go is the non-interacting Green's function, W, S, 
and G' are the grand-canonical potential, self energy, and 
Green's function of the reference system, respectively. In 
this paper, a L = 10 sites cluster (Pig. [T^) is chosen as 
a reference system, and is solved exactly by Lanczos di- 
agonalisation method to obtain the reference self-energy. 
All our calculations are performed at zero temperature 
for the well-known difficulty of including the tempera- 
ture effect into Lanczos method. Since we are looking for 
orbital ordering, a orbital ferromagnetic or antiferromag- 
netic field is used as a variational parameter, in addition 
to the cluster on-site energy. The latter is necessary in 
order to obtain a thermodynamically consistent particle 
density 

The trace in Eq. ^ implicitly contains a sum over 
Matsubara frequencies which needs to be carried out with 
high accuracy. In connection with a Lanczos diagonalisa- 
tion of the cluster Hamiltonian this can be done by means 
of the sum over the single-particle excitation energies ob- 
tained by the band Lanczos [3^ method, as explained in 
Rcf. [35I (see also Ref. [l^). Alternatively, the same accu- 
racy can be obtained more efficiently by an integration 
over an appropriate contour of the complex frequency 
plane, as discussed in Appendix [XI Notice that although 
the contour (see Pig. [6]) mainly runs at a finite distance 
S from the real axis in order to avoid sharp structures 
in the spectral function in the 6^0 limit, the proce- 
dure is exact. There is no need to carry out a 5 — ^ 
extrapolation: this is exactly contained in the additional 
contributions from the "vertical" paths ( C3, C3, C,5, C5, 
C7, Cj in Pig. E]) (see App.E]for details). 



III. RESULTS 



A. Filling dependent single-particle spectral 
function 



In order to gain insight on the physical properties of 
the spinless p-band Hubbard model (Eq. dT])), in this sec- 
tion, we calculate its single-particle spectral function us- 
ing VCA. The VCA has been shown to be an effective 
method to evaluate the single-particle [2^, [sj and two- 
particle nil spectral function of Hubbard-like models. 
The filling-dependent spectral function of the px orbital 
with interaction U = Sty is displayed in Fig. [51 By sym- 
metry reasons, the spectrum of the Py orbital in the non- 
ferromagnetic phase is obtained by interchanging with 
ky. The spectra at different fillings are obtained in the 
respectively stable phase, according to the phase diagram 
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FIG. 2; Single-particle spectral function A(k,ij) of px orbital 
at various fillings From (a) to (c) the fillings of the system are 
1.0, 0.8, and 0.6, respectively. The interaction is U = Sty and 
the transverse hopping is t± — — O.OSty. 



displayed in Fig. |4] (see Sec. IIII B|) . The spectrum of the 
Px orbital is almost fc-independent in the ky direction of 
the Brillioun Zone (BZ), since the dispersion is nearly ID. 
This is quite obviously due to the fact that the transverse 
hopping {t± = — 0.05t|[) is very small. For all fillings, we 
can clearly recognize the upper and lower Hubbard bands 
with a gap of the order of U. 

At half-filling, the spectrum has a ladder structure (see 
Fig. 12^.), which is also characteristic of the t — model 
[2l|, 123, |33] . However, the spectrum is slightly dispersive 
in the kx direction of the BZ, that is, a hole or particle 
is not localized but moves coherently through the lat- 
tice. The small dispersion can be explained by including 
a three-sites term in the t — Hamiltonian [2^ . More 
spectra at half-filling and for different interactions are 
given in Fig. [31 The gap between upper and lower Hub- 
bard bands decreases as the interaction decreases. At 
the same time, the bandwidth become larger because a 
hole (or particle) is easier to move when the interaction 
is smaller. 

Away from half-filling, the quasi-particle spectrum be- 
comes strongly dispersive (see Fig. The shape of the 
spectrum is similar to that of ID free particles, but with 
a strongly renormalized bandwidth. The bandwidth be- 
comes larger and larger when going away from half-filling, 
which means that particles can move easier. Another 
feature that can be seen in Fig. [2] is the spectral weight 
transfer phenomenon between the upper and lower Hub- 
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FIG. 3: Single-particle spectral function A(k,u}) of the px 
orbital at half-filling and for different values of the interaction 
U. Specifically, we have U = 10t|| (a), 6t|| (b), and Sty (c). 



bard bands, which is also observed in the usual single- 
band Hubbard model [l^]. For fillings below half- filling 
(Fig.[2l3,c), the spectrum is transferred to the lower Hub- 
bard band. This is because at low density the parti- 
cles have less chance to doubly occupy the same site and 
therefore have a smaller probability to be in the upper 
Hubbard band Of course, for filling above half- filling 
the situation is reversed due to particlc-holc symmetry. 



B. Antiferromagnetic orbital order 

In this section, wc discuss the "antiferromagnetic" or- 
bital order in this model as a function of filling. At 
half-filling and in the strong-coupling limit U 3> ty, the 
Hamiltonian Eq. (P) can be reduced to a superexchange 
t—J^ model with a positive exchange energy J = 2i|/J7 

[lUilH, H^. Therefore, the Mott state favors a staggered 
("antiferromagnetic") orbital order (see Fig. [Us). To 
study the orbital order within VGA, we add a "virtual" 
staggered orbital field, H'^^p = h'j^p X)r('^r,x — ?^r.i/)e**^ '" 
with Q = (tt, tt), to the cluster Hamiltonian H' . As 
explained in Sec. |TT1 this term is then subtracted pertur- 
batively, and the coefficient is determined by optimizing 
the grand-canonical potential Eq. The correspond- 
ing staggered orbital order parameter, m = X]r(("'r,s) — 
(nr.i,))e**^ '", is then calculated and plotted in Fig. 0^ as 
a function of the interaction U . One can see that the 
order parameter m, which is non-zero for any finite J7, 
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FIG. 4: (Color online) (a), Staggered orbital order parame- 
ter m (□) as a function of interaction U at half-filling, (b), 
Orbital polarization P a.t U — 8f|| (o) and staggered orbital 
order parameter m at f/ = Sty (A) and U = Sty (*) as a 
function of filling n. 

is increasing as U increases and approaches unity in the 
strong coupling limit. This result supports the existence 
of the "antiferromagnetic" orbital order at half-filling. 

The "antiferromagnetic" orbital order is destroyed by 
doping the system away from half-filling. This is illus- 
trated in Fig. Hb, where the staggered orbital order pa- 
rameter m is plotted as a function of filling n at different 
interactions U = 8t|| (denoted by A) and U — 3t|| (de- 
noted by *). Fig. shows that the order parameter 
m decreases sharply when n decreases, and disappears 
completely (m = 0) at fillings n w 0.85 and n w 0.94 for 
U = 8t|| and U ~ 3i||, respectively. We conclude that 
the "antiferromagnetic" orbital order at large U is more 
difficult to destroy than at small U. After m = 0, the 
system enters a featureless "paramagnetic" orbital state. 



C. Orbital order at low filling 

As shown in the spectrum (see Fig. [21 [3]), the spinless 
p-band model has strong ID character for each orbital 
due to the anisotropic hopping, and therefore its band 
structure has a Van Hove singularity near the band edge 
. It is interesting to see whether or not this singularity 
can produce "ferromagnetic" orbital order at low filling 
(shown schematically in Fig. [T}:) [4l|. 

The question is subtle because the same Van Hove sin- 
gularity present in the ID single-band Hubbard model at 
low filling is not sufficient to obtain ferromagnetism. In 
particular, an Hartree-Fock argument provides the wrong 
conclusion that the ferromagnetic state should be lower 
in energy than the paramagnetic state at sufficiently low 
densities and large U in one dimension. This results is 
indeed contradicted by the rigorous Lieb-Mattis theorem 
[i^ . which excludes ferromagnetism for the ID Hubbard 
model, as well as by an accurate analysis based on the 
Gutzwiller wave function. 

While the Lieb-Mattis theorem does not apply to the 
present p-band model, we investigate here whether or not 



an instability of the totally polarized ferromagnetic state 
towards a variational, less polarized, and, ultimately, 
paramagnmetic wave function can be found for the p- 
band model in the low-density and U —^ oo limit. As 
trial wave functions for the less polarized state we use 
the Gutzwiller wave function, as well as a more general 
one, i. e. with lower energy. Despite this, we find that 
the totally polarized state, which, of course, can be solved 
exactly for an onsite interaction, always has the lowest 
energy. While we were not able so far to prove that 
the totally polarized state is the most stable one at suf- 
ficiently low filling, the fact that we have used a quite 
general trial wave function makes us confident that there 
should be no wave functions with a lower energy than the 
totally polarised state. 

We consider a p-band model with N particles in a finite 
L X L square lattice with periodic boundary conditions 
(PBC). For simplicity, we take t± = and U = oo. Quite 
generally, we can expect that if the ferromagnetic phase 
has a lower energy with a finite gap to the paramagnetic 
state for these values of t± and U, its stability region 
should extend to some finite t± and U. 

li N < L, it is quite clear that the lowest energy is 
obtained by putting all particles in the same orbital (say, 
Px) on different "rows". In that case, each particle moves 
independently on its row, so that the kinetice energy is 
minimal and the interaction energy is zero. However, 
this cannot lead to the conclusion that the ferromagnetic 
state is stable at sufficiently low but finite density, since 
for N < L the density vanishes in the thermodynamic 
limit. The crucial question is what happens for N = 
L + 1, i.e. is it more convenient energetically to put the 
next particle in one of the already occupied rows in the 
Px orbital, or to put it in a "column" in the Py orbital? 
If the particle is added to the Px orbital, the system is 
still in a full ferromagnetic state, and the energy change 
AEi of this state with respect to the ground state with 
N = L, \L)x, (which has energy Eq = — 2t||i), is given 
by 

A£;i =-2t||Cos(^)«-2t||(l-^) , (3) 

where, in the last term, we have taken the large-i limit. 
If we add the particle to the py orbital on one of the 
columns (no matter which one), the lowest-energy state 
cannot be determined exactly. Therefore, we approxi- 
mate it by a trial wave function. The simpliest one is the 
Gutzwiller wave function 

1^) - - n,,.n,^y)dl^^^^^^^^ JL)x (4) 

r 

where ^ ^ ^ creates a particle on py orbitals on "col- 
umn" (x) with y wave vector vector qy. The energy in- 
crease can be easily evaluated as 

Ai?2 = -2i||(l--^) (5) 
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Clearly, this energy is larger than ([3]). The reason for 
this is that for the row where the Py particle sits, the 
Gutzwiller wave function of the Px particle has a sharp 
jump at the position of the Py particle, see Eq. This 
leads to an increase of 2t\\/{L — 1) in the kinetic energy 
of the Px particle. A natural improvement consists in re- 
placing the wave function of this row by a smooth sine 
function, sin(7ra;/L), which has much smaller energy in- 
crease. However, the overlap of this sin fun ction with the 
wave function of the other rows, which is yJTjL, is small, 
resulting in a large kinetic energy of the Py particle. A 
proper choice is to linearly combine these two functions. 
This leads to the trial wave function: 



- t 



(6) 



n / 



V2 bra-n ^'^'^Cj^) 



,J0) 



where (^m) explicitly denotes the 2D coordinates of a 
position r, and l,m,n = 0, ,L — 1. Since we are 
using PBC, our assumption that the Py particle is in the 
0th column docs not lead to a loss of generality. The 
coefScients am-n and bm-n can be chosen, for simplicity. 



to be ttm-r 



2 (n-m)j 



and bn 



,^2 {n-m)TT 



This choice docs not affect our conclusions, as discussed 
below. The energy increase for the state in ^ is given 
by: 



2t||(l-- + 0(L-2)) 



(7) 



where a k, 0.5 is a constant. A comparison of the energies 
of this trial wave function with the fully polarized ferro- 
magnetic state, whose energy increase is — 2t|| +0(i^^), 
shows that the latter has a lower energy. Notice that 
more general forms of the coefficients a„i_„ and bm-n do 
not change this conclusion, as they merely modify the 
coefficient a in Eq. ([7]), which, however, remains nonzero 
and positive. 

The above results show that, althought the spinless p- 
band model has strong ID character, it is different from 
the ID Hubbard model. This can be seen by constructing 
a trial wave function in a similar way to Eq. ([6]) for the 
ID Hubbard model with two particles [i^. In this case, 
the total energy for PBC is -4t|| + C/L'^ + 0{L~'^) with 
C = 27r^i|| for the paramagnetic state and C = 47r^i|| 
for the ferromagnetic state [i^, i.e., the ferromagnetic 
state is unstable. This situation is quite different from 
the partly polarized state of the p-band model presented 
above with one particle in the Py orbital and L particles 
in Px orbital. In this case, the motion of the py particle 
is hindered by, and, at the same time, affects the motion 
of the other L particles in the Px orbitals. This leads to 
a much larger energy increase than in the fully polarized 
ferromagnetic state. 

After having discussed the stability of the ferromag- 
netic phase from a more accurate point of view, we return 
to the results of the VC A approximation in the low-filling 



(a) 



(b) 












(c) 



(d) 












FIG. 5: (Color online) Momentum distribution n(fc) for dif- 
ferent values of the interaction U and of the filling n, with 
t± = -0.05*11 . (a), U = 10i||, n = 1.0; (b), U = Sty, n = 1.0; 
(c), U = 10t||, n = 0.6; (d), U = 3f||, n = 0.6. 



region. In the ferromagnetic case it is necessary to intro- 
duce a different on-site energy between the two orbitals as 
a variational parameter. This is equivalent to using the 
cluster chemical potential and a "ferromagnetic" field. In 
the fully polarized case, the saddle point is given by the 
on-site energy of the empty orbital approaching infinity. 
To describe the ferromagnetic phase, we evaluate the or- 
bital polarization P = {nx — ny)/{nx + ny), where Ux and 
Uy are the average occupations of the px and py orbitals. 
Results for P as a function of filling for J7 = 8i|| are plot- 
ted in Fig. [4]d (denoted by o). The orbital polarization P 
is calculated at the respectively stationary point of Q in 
each phase. As in Fig. P vanishes at half-filling and 
remains essentially zero down to a filling of ti sa 0.6. For 
n < 0.6, P rapidly increases as n decreases, and rapidly 
saturates (P = 1) at n « 0.38 indicating a full ferromag- 
netic orbital order state. 

We should stress that one must be careful when in- 
terpreting the VCA results at low filling. First, we can- 
not exclude that finite-size effects, originating from the 
limited size of the reference cluster, could affect the or- 
dered state found in our calculation. This could be the 
case when the exact self-energy is long ranged, so that 
it cannot be accurately described by the self-energy of 
a small reference system. Second, the density obtained 
by VCA shows small discontinuities when the reference 
cluster changes its particle number Therefore, it is 
difficult for VCA to determine the exact critical point for 
the onset of ferromagnetic orbital ordering as a function 
of filling. 

Summarizing this section, our combined VCA and vari- 
ational results are a strong indication, although not a 



6 



proof, for the presence of an orbital ferromagnetic state at 
fow-density and sufficiently large U in the p-band model. 
An exact proof for the absence or existence of ferromag- 
netism at low densities in the p-band model (similarly to 
the Hubbard model [4l|) would be welcome. However, it 
is beyond the goal of the present paper. 

D. Momentum distribution 

In this section, we present results for the momentum 
distribution, as this quantity is directly accessible exper- 
imentally [201, and can be used to detect the possible 
occurrence of orbital ordering. Our results arc shown in 
Fig. m For half-filling and small values of the interaction 
U (Fig. [5Jd), the momentum distribution can approxi- 
mately be seen as the superposition of two ID nonin- 
teracting gases traveling in the two directions x and y. 
Double occupations are present in k space in the middle 
square region of the Brillioun Zone. While double occu- 
pation is allowed for small U, it is strongly suppressed for 
strong interactions. Therefore, for large U , the momen- 
tum distribution is flattened and covers the whole BZ 
(see Fig. [5^). By decreasing the filling away from half 
filling the suppression of double occupation weakens, as 
can be seen in Figs. [SJ:, and d for n = 0.6. The reason is 
that the system has a strongly dispersive spectrum (see 
Fig. [2j;) and, therefore, it is no longer in a Mott state. 
The fact that particles can move quite free in the lattice, 
gives rise to the possibility of double occupation even at 
large interactions. 

Finally, we briefly discuss the experimental signatures. 
The momentum distribution of fermions in the excited 
p-band (see Fig. [S]) is different from that of fermions 
in the lowest s-band [20| . e.g., in the weak interacting 
regime. This can be directly observed in the time of flight 
(TOF) images. The antiferromagnetic orbital order can 
be detected by analyzing the noise correlation function 
from TOF images. In the noise correlation spectrum, the 
s-band fermions produce the antibunching dips at the 
usual reciprocal wave vector of square lattice [i^, H^. 
However, the p-band fermions in the antiferromagnetic 
orbital order state contribute new dips at the reciprocal 
wave vector of doubled unit cell [2l| . 

IV. CONCLUSION 

In summary, we have studied a model for spinlcss p- 
band fermions in optical lattices using the Variational 
Cluster Approach, and, partly, a variational wave func- 
tion. We have computed its single-particle spectral func- 
tion in a wide range of fillings and found a strongly dis- 
persive spectrum at incommensurate fillings. By calcu- 
lating the staggered orbital order parameter, we showed 
that the system is in a staggered ("antiferromagnetic") 
orbital state at half-filling, which is destroyed by dop- 
ing and evolves into a paramagnetic state. In the low- 



density limit and for U = oo, we studied the stability 
of a fully-polarized ferromagnetic state by constructing 
a trial wave function, which extends the Gutzwiller trial 
state. In contrast to the one- and two-dimensional Hub- 
bard model we did not find an instability of the ferro- 
magnetic state towards the paramagnetic solution. In 
particular, for the trial wave function of Eq. ^ (which 
is more general than the Gutzwiller wave function), the 
ferromagnetic state is lower in energy than the param- 
agnetic one. Finally, we have computed, by VGA, the 
momentum distribution and studied its evolution as a 
function of interaction and filling. 
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APPENDIX A: FREQUENCY INTEGRATION 

The sum over Matsubara frequencies in Eq. @ can 
be carried out either (i) dire ctly (see also Ref. [43|) or 
(ii) by the usual procedure [4g| of distorting the con- 
tour to the real axis. For a numerical sum, and es- 
pecially for the corresponding integral at T = 0, both 
procedures present their advantages and disadvantages. 
In case (i) one should first extract the asymptotic (for 
large iw„) part of the integrand and carry out the cor- 
responding sum/integral analytically. In case (ii) there 
is no such problem, as the contribution to the integral 
at an infinitesimal distance (5^0 from the real axis is 
nonzero only within the region where the spectral func- 
tion is nonzero. However, due to the pole structure of 
the integrand, one has to take a finite S for numerical 
purposes. This reduces the precision and introduces ad- 
ditional complications coming from the fact that at large 
Lu the integrand goes like 

The best solution is to distort the integral to the con- 
tour indicated in FiglS) For a sum over Matsubara fre- 
quencies LUn = 27TT{n+ i) of a function g{z) of the com- 
plex variable z, which is analytic everywhere except on 
the real axis , we have 

T ^ e-"«+g(*o;„)=-— i e^°^ fF{z)g{z)dz . 

71 — — OO ^ 

(Al) 

Here, C is the usual contour of the complex plane encir- 
cling the Matsubara frequencies «u;„, /f(-z) = (exp y -|- 
1)~^ is the Fermi function, and 0-1- is a positive in- 
finitesimal. With the usual conditions that g{z) 
for 1^1 OO, and that 

giz*)=9izr, (A2) 
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FIG. 6: Contour in the complex plane in which the frequency 
integration is carried out. 



we can further distort the contour C to the contour in- 
dicated in Fig. [HI Here, Ci and C[ are semicircles at 
infinity, so that their contribution vanish, C2, Cg, and 
6*2, Cg are infinitesimally close to the real axis, as in the 
usual procedure. C5 and C5 are infinitesimally close to 
the imaginary axis, while C4, Cg, and C4, Cg can have 
an arbitrary finite distance 5 from the real axis. We take 
f2 as some upper limit of the spectrum, i. e.. 



Img{uj + iO^) = for any 
By calling 

F{z)^g{z)fF{z) 



> n 



(A3) 



(A4) 



we can first write the contributions to Eq. (jAf P from the 
"horizontal" paths: 

Sh = -- [ ImF{z) dz . (A5) 

JC2+C4+C6+CS 

Notice that for (5 ^ one recovers the usual expression 
(ist , and there is obviously no contribution from the "ver- 
tical" paths. The contributions from C2 and Cg vanish, 
due to Eq. (|A3p (/f(z) is analytic across the real axis). 
Therefore, we are left with 



St, 



1 



Im [g(w -t- i5)fF{oJ + i5)] duj . (A6) 



The advantage of taking a finite d is that the integrand 
is smooth and one only needs few uj points in the numer- 
ical integration in order to achieve a good accuracy, in 



contrast to the conventional case of small 5. For T 
Eq. (|A6p reduces to 



Sh 



1 



Xm g{uj + i6)duj 



(A7) 



The rest of the integral is given by the "vertical" paths, 
for example the contribution from C3, C3 is given by 



1 

'2^ 



1 



F{-D. + ix) i dx- / F{~n ~ ix) {-i) dx 



(A8) 



neF{-n + ix)dx 



and similarly for the contribution from C7, C^; 



1 



Cr+C' 



- I neF{n + ix)dx 

TT 



(A9) 



The latter contributions vanishes for T = or can be 
made exponentially small by taking fl/T ^ 1. The con- 
tribution from the "central" vertical paths C^fi'r^ is sim- 
ply given by the original sum over Matsubara frequen- 
cies, however only for \u}n\ < <5 (we must be wise and 
choose S not to coincide with a Matsubara frequency for 
T > 0). Denoting by ujn^^^ the corresponding maximum 
frequency, we have 



1 



2T ^ TZegiicUn) 



(AlO) 



n=0 



which for T = becomes 



TZe g{ix)dx . 



(All) 



The contributions Eq. (|A8|) . Eq. (jA9|, Eq. (jAlOp arc the 
additional integrals to be carried out to compensate for 
the nonvanishing value of 5. We stress that the result is 
exact for any (even large) value of 6 > 0. The numeri- 
cal advantage is that the integrand is everywhere smooth 
except for small temperatures and on the path C5 near 
u!n = whenever g{uj) has poles close to a; = 0. More- 
over, all integrals are carried out in a finite domain, so 
there is no need to carry out extrapolations. 
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